rm(list = ls())

library(readr)
library(ggplot2)
library(stargazer)
library(lme4)
library(dplyr)
library(haven)
library(texreg)
library(tidyr)

load("Final_JPR/data/paperdata.RData")
data <- paperdata %>%
       filter(tek_count>0)

source("Final_JPR/Code/coefplot.R")
source("Final_JPR/Code/marginal_effects.R")
source("Final_JPR/Code/setup.R")


######### Figure 4


data <- data %>% dplyr::mutate(peaceyears_terr_sq = peaceyears_terr^2,
                               peaceyears_terr_cub = peaceyears_terr^3,
                               peaceyears_gov_sq = peaceyears_gov^2,
                               peaceyears_gov_cub = peaceyears_gov^3)

## dependent variable
dv <- 'onset_do_terr_flag' 
## controls
ctrs <- c('lineq_total', 'status_excl','lsize','family_warhistdummy','family_downgraded2', 
          'ctygdppc_log', 'ctypop_log', 'peaceyears_terr','peaceyears_terr_sq','peaceyears_terr_cub')

ivs_1 <- c('best_tlineq_total', 'tekbest_status_excl')
ivs_2 <- c('worst_tlineq_total', 'tekworst_status_excl')
ivs_3 <- c("median_tlineq_total",'tekmedian_status_excl')
ivs_4 <- c('nearst_tlineq_total','teknearst_status_excl')


f_do_1 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_1, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_2 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_2, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_3 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_3, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_4 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_4, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))

#
ivs <- c('best_tlineq_total','worst_tlineq_total','median_tlineq_total','nearst_tlineq_total')


## using log of transnational inequality
model_onset_terr<- list()
for (i in 1:length(ivs)){
       model_onset_terr[[i]] <- glm(eval(parse(text = paste0('f_do_', i))), 
                                    data = data, family = binomial(link = "logit"))
}



dv <- 'onset_do_gov_flag' 


ctrs <- c('lineq_total', 'status_excl','lsize','family_warhistdummy','family_downgraded2', 
          'ctygdppc_log', 'ctypop_log', 'peaceyears_gov','peaceyears_gov_sq','peaceyears_gov_cub')

ivs_1 <- c('best_tlineq_total', 'tekbest_status_excl')
ivs_2 <- c('worst_tlineq_total', 'tekworst_status_excl')
ivs_3 <- c("median_tlineq_total",'tekmedian_status_excl')
ivs_4 <- c('nearst_tlineq_total','teknearst_status_excl')


f_do_1 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_1, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_2 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_2, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_3 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_3, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_4 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_4, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))

#
ivs <- c('best_tlineq_total','worst_tlineq_total','median_tlineq_total','nearst_tlineq_total')


## using log of transnational inequality
model_onset_gov<- list()
for (i in 1:length(ivs)){
       model_onset_gov[[i]] <- glm(eval(parse(text = paste0('f_do_', i))), 
                                   data = data, family = binomial(link = "logit"))
}


## dependent variable
dv <- 'incidence_terr_flag' 
## controls
ctrs <- c('lineq_total', 'status_excl','lsize','family_warhistdummy','family_downgraded2', 
          'ctygdppc_log', 'ctypop_log',  'peaceyears_terr','peaceyears_terr_sq','peaceyears_terr_cub')

ivs_1 <- c('best_tlineq_total', 'tekbest_status_excl')
ivs_2 <- c('worst_tlineq_total', 'tekworst_status_excl')
ivs_3 <- c("median_tlineq_total",'tekmedian_status_excl')
ivs_4 <- c('nearst_tlineq_total','teknearst_status_excl')

f_do_1 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_1, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_2 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_2, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_3 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_3, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_4 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_4, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))

#
ivs <- c('best_tlineq_total','worst_tlineq_total','median_tlineq_total','nearst_tlineq_total')


## using log of transnational inequality
model_incidence_terr<- list()
for (i in 1:length(ivs)){
       model_incidence_terr[[i]] <- glm(eval(parse(text = paste0('f_do_', i))), 
                                        data = data, family = binomial(link = "logit"))
}



dv <- 'incidence_gov_flag' 


ctrs <- c('lineq_total', 'status_excl','lsize','family_warhistdummy','family_downgraded2', 
          'ctygdppc_log', 'ctypop_log', 'peaceyears_gov','peaceyears_gov_sq','peaceyears_gov_cub')

ivs_1 <- c('best_tlineq_total', 'tekbest_status_excl')
ivs_2 <- c('worst_tlineq_total', 'tekworst_status_excl')
ivs_3 <- c("median_tlineq_total",'tekmedian_status_excl')
ivs_4 <- c('nearst_tlineq_total','teknearst_status_excl')

f_do_1 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_1, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_2 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_2, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_3 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_3, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_4 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_4, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))

#
ivs <- c('best_tlineq_total','worst_tlineq_total','median_tlineq_total','nearst_tlineq_total')


## using log of transnational inequality
model_incidence_gov<- list()
for (i in 1:length(ivs)){
       model_incidence_gov[[i]] <- glm(eval(parse(text = paste0('f_do_', i))), 
                                       data = data, family = binomial(link = "logit"))
}


#model_onset_terr, model_onset_gov,model_incidence_terr,model_incidence_gov

p1 <- marginaleffect_logit2(ModelResults = model_onset_terr, n.sim = 1000, 
                            varname = c("best_tlineq_total", "worst_tlineq_total",
                                        "median_tlineq_total","nearst_tlineq_total"),
                            data = data,
                            dv = 'onset_do_terr_flag',
                            val1 = 0, val2 =1, clusterid = "gwgroupid") 


p2 <- marginaleffect_logit2(ModelResults = model_incidence_terr, n.sim = 1000, 
                            varname = c("best_tlineq_total", "worst_tlineq_total",
                                        "median_tlineq_total","nearst_tlineq_total"),
                            data = data,
                            dv = 'incidence_terr_flag',
                            val1 = 0, val2 =1, clusterid = "gwgroupid") 


p3 <- marginaleffect_logit2(ModelResults = model_onset_gov, n.sim = 1000, 
                            varname = c("best_tlineq_total", "worst_tlineq_total",
                                        "median_tlineq_total","nearst_tlineq_total"),
                            data = data,
                            dv = 'onset_do_gov_flag',
                            val1 = 0, val2 =1, clusterid = "gwgroupid") 


p4 <- marginaleffect_logit2(ModelResults =model_incidence_gov, n.sim = 1000, 
                            varname = c("best_tlineq_total", "worst_tlineq_total",
                                        "median_tlineq_total","nearst_tlineq_total"),
                            data = data,
                            dv = 'incidence_gov_flag',
                            val1 = 0, val2 =1, clusterid = "gwgroupid") 




df_logit <- bind_rows(p1, p2, p3, p4)


p = ggplot(df_logit , aes(colour = dv)) + 
       geom_hline(yintercept = 0, lty = 2) +
       geom_linerange(aes(x = id, ymin = low,
                          ymax = high),
                      lwd = 1, position = position_dodge(width = 1/2)) +
       geom_pointrange(aes(x = id, y =median, ymin =low,
                           ymax =high, shape = dv),
                       fatten = 7,
                       lwd = 1, position = position_dodge(width = 1/2),
                       fill = "WHITE")+
       geom_text(aes(x = id +0.15 ,
                     y = median,
                     label = format(round(median, 2))),position = position_dodge(width = 1/2)) +
       theme_bw() + xlab("") + ylab("First Difference in Predicted Prob.") + 
       theme(legend.position="bottom",
             legend.title=element_blank(),
             legend.text = element_text(colour="blue", size=10),#black
             legend.key.height = unit(1, 'cm'),                     
             axis.text= element_text(size=11),
             text = element_text(size=10),
             plot.title = element_text(hjust = .5, size = 12, face = "bold"))+
       ggtitle("")+
       #scale_colour_manual(values = c("black", "black","black", "black"),labels = rev(c("Conflict Onset (territorial)", "Conflict Onset (governmental)",
        #                                                               "Conflict Incidence (territorial)","Conflict Incidence (governmental)")))+
       scale_colour_discrete(labels = rev(c("Conflict Onset (territorial)", "Conflict Onset (governmental)",
                                            "Conflict Incidence (territorial)","Conflict Incidence (governmental)")))+
       scale_shape_manual(values = c(13, 15,16,17),labels = rev(c("Conflict Onset (territorial)", "Conflict Onset (governmental)",
                                                                "Conflict Incidence (territorial)","Conflict Incidence (governmental)")) )+

       guides(color = guide_legend(nrow = 2),
              shape = guide_legend(nrow = 2, byrow = TRUE,
                                   override.aes = list(size = 1))
       )+
       
       scale_x_continuous(breaks = 1:4,labels =parse(text =rev(c( 
              "Nearest~TEK",
              "Median~TEK",
              "Worst~TEK",
              "Best~TEK")))) 

ggsave("Final_JPR/Figures/fig4.pdf", width = 7, height = 4.5, units='in', dpi=600)
ggsave("Final_JPR/Figures/fig4.jpg", width = 7, height = 4.5, units='in', dpi=600)
ggsave("Final_JPR/Figures/black_figures/fig4.jpg", width = 7, height = 4.5, units='in', dpi=600)

